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§ ■ ABSTRACT 
(N 

<D • We examine whether massive binary black holes in spherically symmetric 



bulges of galaxies can achieve coalescence through the emission of gravitational 



^ ■ radiation under the action of stellar dynamical processes alone. In particular, 

we address the importance of the Brownian motion of a binary's center-of-mass 

^ ■ to its continued interaction with stellar orbits that allows it to keep hardening: 

^ ■ the restoring force holding the binary at the center decreases as it depletes the 

lO ■ central stellar density, causing the amphtude of wandering to increase. We use 

^ ■ an analytical model and N-body simulations to calculate the time required to 

. reach the gravitational-radiation dominated stage. We find that a substantial 

. fraction of all massive binaries in galaxies can coalesce within a Hubble time. 

o 

^: 

! celestial mechanics — galaxies: nuclei — methods: N-body simulations 



Subject headings: black hole physics — galaxies: kinematics and dynamics 



. 1. Introduction 

Observations indicate that most galaxies harbor supermassive black holes (BHs) at 
their centers (e.g., Magorrian et al. 1998). Since galaxy mergers are inevitable, there arises 
the question of the dynamical fate of the binary black holes (BBHs) that result from these 
mergers. This problem was first addressed by Begelman, Blandford & Rees (1980) (for 
recent work, see, e.g., Milosavljevic & Merritt 2001, 2002, Yu 2002). Dynamical friction 
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causes the black holes in the original galaxies to sink to the center of the merged galaxy, 
and then become bound in a binary system. The semi-major axis of the binary, a, continues 
to shrink owing to dynamical friction until it reaches the point at which the binary becomes 
hard: ah — Gm2/4o"^, where a is the one-dimensional stellar velocity dispersion, G is the 
gravitational constant, and m2 (< mi) is the mass of the lighter of the two black holes 
(following Quinlan 1996). After this stage, the dominant mechanism by which the binary 
shrinks and loses energy is by three-body interactions with stars that pass within a distance 
~ a of the center of mass of the binary. The binary transfers energy to the stars by ejecting 
them at a much higher velocity. Once it shrinks to the point at which energy loss to 
gravitational waves becomes dominant, the binary can quickly coalesce. 

Begelman, Blandford & Rees (1980) conjectured that hardening of a tight binary 
by three-body interactions would continue only until it had interacted with all the stars 
contained within its loss cone (Prank & Rees 1976); if this was not sufficient to drive it to 
the gravitational wave emission stage, the binary would stall at a certain separation, and 
harden very slowly at the two-body relaxation timescale (the timescale at which the loss 
cone can be repopulated), which turns out for many actual galaxies to be longer than the 
Hubble time. The lack of observational evidence for the existence of massive BBHs has led 
to suggestions that BBH merging may be facilitated by other processes, such as the loss 
of angular momentum to gas (e.g., Begelman, Blandford & Rees 1980, Gould & Rix 2000, 
Armitage & Natarajan 2002), or triaxiality of the host galaxy (Yu 2002). 

It is therefore of interest to inquire whether stellar dynamical processes alone would be 
sufficient to drive the BBH to the gravitational wave emission stage. Quinlan & Hernquist 
(1997) pointed out the importance of the wandering of the binary to its hardening. In this 
paper, we investigate the shrinking of a hard black hole binary in a spherically symmetric 
galaxy, and ask particularly how Brownian motion of its center of mass can contribute to its 
continued hardening beyond the stage of loss cone evacuation. In §2, we set up an analytic 
model that appears to account for the results of N-body simulations described in §3. We 
comment on the impact of BBH wandering in §4, and derive characteristics of the merger 
in §5, on the assumption that hardening proceeds at the rate deduced in §2 and §3. We 
summarize our findings in §6. 

2. The Model 

Consider two black holes of masses mi and 7712- Let fiu = mim2/{mi + 7722) be the 
reduced mass of this BBH system. Also, let m2 be the lighter of the two BHs, with the 
mass ratio being q = 1712 /mi < 1. 
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For illustrative purposes, we take the host cluster of stars to be described by a Plummer 
model of total mass M and length parameter tq. Thus, the density and potential profiles 
are given, respectively, by 

. ^ _ 3Mrg 1 
W 4^ (^2 + ^2)5/2' 

where G is the gravitational constant and r is the radial distance from the center of the 
stellar system, which is taken as the origin. In this case, the one-dimensional stellar velocity 
distribution in the core is given by cr^ = GM/6ro. 

The phase space distribution function of the stellar system, /, depends in general both 
on position r, and stellar velocity iTa, and is defined such that f{f,iTa) d^fd^Va is the mass 
in stars in the phase space volume d^fd^Va- For the spherically symmetric Plummer model, 
we assume / to be a function of the relative energy per unit mass £ only (and independent 
of specific angular momentum), where £ — —\v1 — $(r) = ^(r) — \v1, ^{r) — — $(r) being 
the relative potential. The distribution function is then (e.g., Binney & Tremaine 1987) 

Consider an encounter between the BBH and a star of mass m*. Let the impact 
parameter be p and the initial relative velocity (when the star is at infinity) be Vq. If 
Tmin is the minimum distance reached by the star to the BBH, and Vmax is the velocity at 
that point, then the angular momentum is / = pVo = Tminymax- The reduced mass of the 
BBH-star system is ^ ~ m*, assuming the BHs to be much more massive than the stars. 
Then by conservation of energy. 



1 2 G{mi + m2)m^ 1 2 
E = -m*VC_ ^ = -m^Vn 

Eliminating Vmax and rearranging. 



iy -*^max ~ r)'"'*''0 (4) 



2_ 2 I 2G{m, + m2) 

y ~ ' min ^ T/2 V^/ 

We set Tmin = a, the semi-major axis of the BBH, implying that an encounter occurs (i.e., 
according to our picture, the star interacts with the BBH and gets thrown out by the 
slingshot effect, while a shrinks) if < + m^^pnsla. 
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The number of stars per unit volume with velocity Va is f{va)(Pva/m^. Hence, the 
number of encounters per unit time is 

lMd\^^p2y^, (6) 

where Vq = |u — Wal, v being the velocity of the BBH. In order to take into account the 
wandering of the BBH's center of mass, this expression must be averaged over v (for a 
similar calculation, see Binney & Tremaine [1987], §8.5). It was shown in Chatterjee, 
Hernquist & Loeb (2002a, b) that the distribution of each component of v is Gaussian, so 
that the distribution of the magnitude, v, is W{v) = T^s-e-'^'/^'^bs with ct^^ = 9rS^+^- 
Thus, averaging over W, 



ten \ph:m 
or 

1 27r 



ten \p2^m. 



where A — y^^3 (oSj^ • Note that in the core, obtains an approximately constant value, 
il) ~ GM/tq. 

The above integral can be simplified using the approximation jiT — t}^! ~ since 
the BBH moves slowly relative to the stars when its mass is much larger than that of an 
individual star. The integration then gives, 

1 327r2AV'^^/V llG(mi + m2)\ 

— — — t;^ + ^ ; (9) 

ten 99m, \ 2 I ^ ' 

The BBH's binding energy is £■ = If each star on average extracts energy A£^, 

where Ai? = G'//i2m*/a, is roughly the potential energy of the BBH-star system at the 
point of closest approach (giving ISEjE — ^^^^ ) , then 

A£; _ Gmim2m^ 32%^ Aip^^^^ f 11 G{mx + m2)\ _ Gmim2 da 

'At ~ mi + m2 99m* ^ T )~ 2a^ It' ^ ' 

and simplifying, 

1 da eAn^A^jj^^/^ ^^^^ 



a2 



(a + iiG(y»2)^ dt 99(mi + m2) ' 



Thus, there are two limiting regimes for the temporal evolution of a: 
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• If a » llG(mi + 1112) /2t/j, then = mi+m^ ' '^^^^^ C is a constant, and we get 
l/aociV2. 

• If a <S llG(mi + m2)/2ip, then ^(^) = a constant, and we get 1/a oc t. 

For intermediate values, we get a power law, l/acxf^, with 0.5 < d < 1.0. 

Now, llG'(mi + m2)/2tl) = ^^^^^^o- The BBH becomes hard at an - f^ro, where 
1712 is the hghter of the two BHs. Thus, after the BBH becomes hard, we always tend to the 
second case above, 1/a oc 

In our numerical experiments (to be described in the next section), we have found 
temporal power-law indices for 1/a which are always bounded between ~ 0.5 (for heavy 
BBHs) and 1 (for light BBHs). The above simphfied theory predicts a power-law index of 
1 for hard binaries. However, our experiments showed that heavy BBHs (which contain 
more than a few percent of the mass of the Plummer galaxy) cause the central density 
of the galaxy to drop quickly, and so it is not appropriate to assume, as we did above, 
that the stellar distribution function remains unaffected by the BBH. In particular, the 
left-hand-side of equation (11) needs to be ^2(~_^ei) '^here Ci is now some undetermined 
term. If Ci is not small compared to a, then it is possible to get 1/a oc t'^, with d < 1. This 
description is suggestive for why we obtained such limiting behaviors, and it is due to the 
form of the expression (5) for the impact parameter being a sum of two terms. 

For low mass BBHs, the distribution function is not much affected by the hardening of 
the binary, and so we may proceed to derive the temporal slope of the 1/a line, 

dt\a) 2l7r V y 
Note that this result is independent of the masses of the two BHs. 

We may compare the above value with the hardening rate presented in other work in 
a slightly different form: ^(^) = where H is & dimensionless hardening rate, po is 

the central stellar density, and a is the one-dimensional stellar velocity dispersion in the 
core. Comparing with the above equation, we obtain H — 9.4. This is very similar to the 
hardening rate observed in simulations by Milosavljevic & Merritt (2001). It also compares 
well with the value if ~ 15 (independent of BH mass ratio and orbital eccentricity), which 
is the hardening rate derived from three-body scattering experiments in the limit of a very 
hard binary (Quinlan 1996, Hills 1983). 

In our simulations, we made the choice Tq = Sir/ 16 and GM — 1, for which one gets 
s — ^(^) = 20.6. Note, however, the approximations that have been made in the derivation. 
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We have assumed that the distribution function of the stars remains unchanged, which is 
clearly false: the central density of the stellar system falls, as shown in the next section. 
For low mass BBHs (say m < 0.01), however, the decline is not as much as for heavy BBHs, 
and we expect the inverse semi-major axis of the BBH to increase at a constant rate, with 
s given roughly by the expression above. 

We have assumed AE/E = 2m^,/(mi + w^), a form which is in reasonably good 
agreement with results of scattering experiments at several mass ratios (see Roos 1988, Hills 
1983). It is, of course, an expression for the average energy extracted by a passing star, and 
need not accurately reflect the actual process of energy extraction for all impact parameters. 
Given that the precise form of this expression suitable for use in the massive BBH case 
is uncertain (for example. Hills & FuUerton [1980] find that for equal-mass binaries, the 
right-hand side of the above expression may be multiplied by 0.7), it is perhaps optimistic 
to expect this derivation to give the correct slope s. 

3. Results from Numerical Experiments 

The numerical simulations were carried out using the SCFBDY program, which is 
described in detail in Quinlan & Hernquist (1997) and Chatterjee, Hernquist & Loeb 
(2002a). In our code, each BH interacts with the other BH and the stars through the 1/r^ 
force, but stars interact with other stars through a mean field expansion of the gravitational 
potential that is updated self-consistently with time (note that the simulation thus does not 
treat star-star relaxation accurately [see Hernquist & Barnes 1990]; this is justified since the 
two-body relaxation timescale in real galaxies is much longer than the Hubble time). The 
particles are moved with individual step-sizes using the fourth-order integrators of Aarseth 
(1994): NBODYl for the stars and NBODYl or NB0DY2 for the BHs. 

The mean field expansion of the gravitational potential is carried out in a set of basis 
functions that includes both spherical and non-spherical terms (Hernquist & Ostriker 1992). 
However, since the galaxy remained nearly spherical during the integration, it made not 
much difference to the hardening of the BBH whether a few or zero non-spherical terms 
were included in the potential expansion. 

We verified that numerical relaxation was not an important factor in our simulations 
by forcing the stars to interact only with a fixed background potential, and not with each 
other, and finding our results to be similar to those obtained using the full SCFBDY 
code. The insignificance of numerical relaxation in our simulations is further supported by 
experiments in which we forced the center of mass of the BBH to remain fixed, as described 
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below. 

Our results are presented in a system of units in which G = M — 1, and Tq = 37r/16, so 
that the energy of the initial galaxy without the BHs is -1/4 (Heggie & Mathieu 1986). The 
stars (of equal mass m* — 1/N, where N is the number of stars) are initially distributed 
according to the Plummer distribution function, and the black holes are started out on 

nearly circular orbits at r and — r, with r = 0.3. (Wc tried a range of initial conditions for 
the BHs, including elliptical orbits; it made no difference to the hardening of the binary.) 
We have done a range of experiments, with binary mass ranging from rui = m2 = 0.00125 
to mi = 1712 = 0.01, and with particle numbers of up to iV = 400, 000. 

Figure 1 shows the results of such an integration with mi — 1712 — 0.00125, and 400,000 
stars. Shown are the separation between the BHs, R, the inverse semi-major axis of the 
BBH orbit, 1/a, the eccentricity of the orbit, e, and the central density, p, of the stellar 
system averaged over a sphere containing a mass in stars somewhat larger than the mass of 
the BBH. The BHs become bound at t = 77.3. The BBH becomes hard once a shrinks to 
the value ah — G'm2/4(7^ = 1.104 x 10~^, or when 1/a/j = 905. As can be seen, the rate of 
change of 1/a is nearly constant after the binary has become hard; namely, 1/a oaf^ with 
d~ 1. 

Figure 2 shows the results of a similar integration with rrii — 1712 — 0.02, and 200,000 
stars. The BHs become bound a,t t — 4.4, and the BBH becomes hard when a — ah — 0.018, 
or when 1/a — 1/ ah — 56.6. The rate of hardening is clearly not constant after this stage. 
Indeed, this is at the opposite extreme from the previous case, with the hardening being 
well fitted by 1/a oct^ with d ~ 0.5. 

For intermediate values of mi and 1712, the hardening occurs at the rate 1/a oc 
with 0.5 ^ d ^ 1.0. We have not found values of d significantly outside this range. This is 
consistent with the explanation given in §2. Notice that for mi = 1712 — 0.00125, the density 
at the position of the BBH center of mass drops gently over the course of its hardening. 
This is expected since the BBH expels stars from the center as it hardens. The magnitude 
of the drop is small, however, and thus to a first approximation, the stellar distribution 
function remains unaffected. On the other hand, the heavy BBH changes the distribution 
function to a greater extent, as is clear from the bottom right panel of Figure 2. Thus, in 
the former case, the theory developed in the previous section holds, and the binary hardens 
at a constant rate; in the latter case, we get a qualitatively different behavior, for the 
reasons outlined in §2. 

The rate of hardening over time is close to constant in three cases for which we ran 
simulations: mi = m.2 — 0.00125, mi = m2 = 0.0025, and mi = m.2 — 0.005. Like Quinlan 
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& Hernquist (1997), we find that the slope of the 1/a fine, given by s — ^(1/a), varies 
with the number of stars. As the number of stars is increased, s falls systematically, until 
there are roughly 200,000 - 400,000 stars, when it stabilizes to a particular value, Sq. This 
decrease in the slope as the number of stars increases has also been observed by Makino 
(1997) and Hemsendorf, Sigurdsson & Spurzem (2002), though not by Milosavljevic & 
Merritt (2001). Quinlan & Hernquist (1997) explained it as resulting from two competing 
factors governing the wandering of the BBH's center of mass: on the one hand, as the 
number of particles, N, rises, the binary wanders less from the center, and thus it reduces 
the density at the center more, and the hardening rate decreases; on the other hand, the 
restoring force of the stellar potential goes down as the central stellar density drops, causing 
an increase in the BBH's wandering. At a particular value of N, there is a balance between 
these two factors, and the hardening rate stops changing for higher values of N. 

Quite apart from the variation of s with N, we found that the limiting value of s 
changes also with the mass of the BBH. We found that for mi — m2 — 0.005, s stabilizes 
to a value of sq ~ 8; for mi = m2 = 0.0025, Sq ~ 11; and for mi = m2 — 0.00125, sq ~ 15. 
Recall that s = 20.6 according to the simplified theory of §2. That value is expected to 
be obtained in the limit that the stellar distribution function remains unchanged by the 
infiuence of the binary. Our numerical experiments indeed reveal that this becomes a better 
approximation as the binary mass is reduced. However, the question remains whether this 
value is actually reached for still less massive BBHs. 

The orbital eccentricity of the BBH in its hard state varies according to the mass of 
the binary. For massive binaries, such as for the case mi = m2 = 0.02 (Figure 2), the BBH 
orbits are circular to a good approximation. For smaller masses of the BBH, the eccentricity 
is larger: e.g., for mi — m2 — 0.00125 (Figure 1), the eccentricity in the hard state is 
roughly 0.3, and trending very slightly upward. These values are consistent with Quinlan 
& Hernquist (1997) and Milosavljevic & Merritt (2001), and with the results of scattering 
experiments in the restricted three-body approximation (Quinlan 1996), according to which 
the BBH eccentricity hardly grows if the initial eccentricity is moderate (cq < 0.3), and does 
not grow by much if it is larger, before the gravitational radiation stage sets in. However, 
Hemsendorf, Sigurdsson & Spurzem (2002) and Aarseth (2002) have recently found that the 
eccentricity may grow to values close to 1, resulting in rapid coalescence by the emission of 
gravitational radiation (see §5). The reason for this discrepancy is not clear, and may have 
to do with the initial conditions chosen for the simulations. 

It is of interest to know the extent of wandering of the BBH. According to the formulae 
developed in Chatterjee, Hernquist & Loeb (2002 a, b), the mean squared position and 
velocity of the center of mass of the BBH in the hard stage should, for the case of the 
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Plummer potential, be given, respectively, by 



(13) 



cm 



3(mi + 777,2) ' 




2GMm^ 



(14) 



3ro(r?7i + m.2) 



This assumes that the stellar distribution function is unchanged, which is not the case for 
a massive BBH; in fact, the wandering should be larger than indicated by these values, 
since the restoring force on the BBH drops as the central stellar density is diluted by the 
BBH. However, the above expressions provide a reasonable description for the case of light 
BBHs. Indeed, we find from simulations that for such cases, the wandering is consistently 
larger than the values provided by the above formulae, but not by more than a factor of 
~ 2 (see Table 1). The reasons for the enhanced wandering could still be the somewhat 
reduced stellar density at the center of the cluster, as well as the enhanced random motion 
of a binary compared with that of a point mass, as discussed by Merritt (2001). Thus, 
the binary is able to interact with a larger region of stellar phase space than calculated 
according to the simple model of §2. 



Quinlan & Hernquist (1997) reported an experiment in which the center of mass of the 
binary was held fixed at the origin by the application of an artificial constraint force: they 
found that hardening halted soon after the binary became hard, following the evacuation 
of the loss cone. We have confirmed this result in our own experiments (see Figure 2 for 
an example): the BBH becomes hard when its semimajor axis reaches Uh, and continues 
hardening by three-body interactions with the stars until a = aic (o/c < ah), when the loss 
cone is empty. These experiments show the importance of wandering in keeping the BBH 
hardening. (Incidentally, the fact that hardening stops in this case is another indication that 
numerical relaxation is not causing spurious refilling of the loss cone in our simulations.) 

According to traditional theory (e.g., Begelman, Blandford and Rees [1980]), hardening 
of the binary should stop after the BBH has interacted with all the stars in its loss cone, 
as long as its wandering is unimportant. However, there is no sign in our experiments that 
a "hole" develops in the density distribution around the binary as it hardens - instead, 
the central density drops slowly as the binary hardens (see Figure 1). On the contrary, as 
explained in the previous section, the wandering of the BBH increases as the central density 
falls, allowing it to interact with new regions of stellar phase space, and to keep hardening 
at a rate consistent with 1/a oc f^. The value of d 1 is reahzed in the hmit where the 



4. The Impact of Wandering on the Hardening of the BBH 
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BBH makes up a small fraction of the host stellar system mass, as is typically the case for 
galactic bulges. 

What if the mean wandering radius of the BBH were smaller than the semi-major axis 
at which the binary becomes hard (a/j)? Then the binary would wander "within itself", 
so to speak; in this case one may wonder whether it would not encounter fresh regions of 
phase space, and perhaps hardening would stop (the bottleneck, if there is one, would occur 
at the stage at which the binary becomes hard, because wandering only becomes more 
important as the semi- major axis of the binary shrinks). Is there a sign of this happening 
in the simulations? 

Unfortunately, to do this experiment in the regime in which wc arc most interested 
(i.e., for BBHs of low mass, so that the hardening rate is constant over time), one would 
require the particle number to be prohibitively high. However, it is possible to perform 
this experiment for more massive BBHs. For example. Table 2 lists the experimental 
root- mean-square value of the binary's Brownian motion, rem, for a number of experiments, 
with the corresponding value of ah', in each of the cases, rem < O'h- Hardening does not stop 
at or near this point, but continues indefinitely beyond this value at the rate 1/a ocf^, with 
0.5 ^d^ 1.0. Figure 2 shows the hardening of a BBH with mi = = 0.02, with 200,000 
particles. The rms value of the binary's Brownian motion is 0.013; the binary becomes hard 
when a — ah — 3m2a/2M = 0.018, or when 1/ah — 56. 

In each of these cases, however, we find that rem > die', moreover, the value of rem 
is very different from the value predicted by equation (13). This indicates that when the 
binary is allowed to wander, hardening no longer stalls at the semimajor axis aic] rather, 
the amplitude of its wandering increases as the central density of the star cluster declines 
(see Figure 3 for an explicit example of this), resulting in a lower restoring force holding the 
binary at the center. This, in turn, allows the binary to keep hardening by interacting with 
new stars that it would not have encountered if it did not wander. Thus, the "loss cone" is 
repopulated at a rate higher than the two-body relaxation rate by virtue of the increased 
wandering of the BBH's center of mass. 



5. Implications for the Coalescence of Massive Black Hole Binary Systems 

We have evidence, therefore, that the BBH can continue to harden beyond the point 
at which the loss cone would be evacuated were the binary's center of mass to be fixed 
at the origin, and that, moreover, it does so at roughly the rate calculated in §2. If we 
assume that hardening can continue until the point at which the emission of gravitational 
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radiation becomes the dominant mode of energy loss for the binary, we can calculate when 
that should occur. 

Prom equation (12), the hardening timescale of the BBH is 



a 



a 



21n /_rl V^h 



T =i:^A7^\ - ■ (15) 



256^2 \GM I a 



The hardening timescale for gravitational radiation is (Peters 1964) 

g _ 5 c'a^ (1-6^)7/2 

~ g ~ 64 6-3^12(7711 + 7712)2 1 + 7362/24 + 3767 96' ^ ' 

where c is the velocity of hght, and e the eccentricity of the BBH's orbit. The dependence 
of the gravitational radiation timescale on eccentricity is weak unless e is close to 1. As 
noted in §3, e is very small for heavy BBHs, but rises to larger values (~ 0.2) for lighter 
BBHs. In the absence of a precise characterization of the variation of e with BBH mass, we 
take 6 = in the above expression; thus, our result for the time of onset of gravitational 
wave emission and coalescence should be regarded as an upper limit, as a higher eccentricity 
would allow this stage to be reached more quickly. 

Equating the two expressions, we have that the semi-major axis at which gravitational 
radiation losses start to dominate, Ggr, is given by 

5 2l7r ^5/2^3^(1 ^ ^)^5/2 

""'"^ 20 V2 c5MV2 ' ^'^> 

where q is the mass ratio q — 7772/7771 < 1. 

From equation (12), assuming that ttgr ^ ah, we can derive the time between the 
binary's becoming hard at semi-major axis g/i, and the time at which 

256^2 GV2MV2a^/ 



or 



256v^V2l7r; GMV5ml^^[q{l + q)]^^' 
The dependences on M, 7771 and q are rather weak. 

Wc can also calculate the total number of stars, Nint, the BBH must interact with from 
the time it becomes hard to the time tgr- We have AE — G//i2777*/g = ^^j^- Therefore 

dE = ^^^dN, or H = ^TTg. But E = G'777 i7772/2g, and so g = ^S^^f- Equating 
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the above, we get ^ = — Integrating, the total number of stars the BBH interacts 

' ° da m» 2o ° °' 

with is 

where a/j = 3m2ro/2M = 3migro/2M. The total mass with which the BBH interacts is 
then simply 



Mint = 



(777,1 + 777,2) (an 



l„(^). (21) 



As will be seen from Figure 4c, the total mass in stars the BBH must interact with is 
roughly of the order of itself. Therefore, according to this theory, the BBH should not 
change the stellar distribution function much, as long as the mass of the BBH is not 
significant in comparison with the total mass in stars. This result is consistent with our 
simulations for small mass BBHs. 

Some of the above relations are depicted for physically reasonable stellar systems 
and BBHs in Figure 4. Instead of presenting results in terms of M and ro, we use more 
physically and observationally relevant variables, such as the one-dimensional stellar velocity 
dispersion, o", and the central density of the stellar system, po- The relations amongst these 
quantities for the Plummer model are as follows: 

3M 



47rro ' 
GM 



6ro 

We then obtain 

5 _ 63V3 GXg(i + gy ...X 

80 c% ' ^ ^ 

and 

- 1024 V63v^J GV^mT[q{l + q)Y/'pT ' 

Figures 4a and 4b show that for physically reasonable masses of the BBH and 
characteristics of the host galaxy, the semi-major axis can shrink to the value a^^ within a 
Hubble time, after which coalescence occurs quickly by the emission of gravitational waves. 
The actual values of agr are very small (Figure 4d). 

The above calculations assume that the binary hardens according to 1/a oc t, which is 
the behavior we see in the simulations when the mass of the BBH is less than about 1% of 
the total mass in stars. Is this a reliable guide to what happens in actual galaxies? 
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The ratio of the mass of the BH to the mass of the galaxy is observed to be quite low 
(on the order of 10~^: see Merritt & Fcrrarcsc 2001, Wandel 1999, Magorrian et al. 1998), 
and therefore the linear hardening behavior is most likely to apply in the cases of practical 
interest. 

However, it is interesting to ask what would happen if the BBH hardened more 
slowly than at a constant rate. As an example, consider the case shown in Figure 2, with 
mi = 7712 = 0.02, and N = 200, 000 stars. The hardening behavior is well fit by the 
expression 1/a = 21.7 x t^'^. For purposes of illustration, let us use a^^.^ to denote the 
semimajor axis at which gravitational wave radiation becomes dominant if we assume that 
the BBH hardens at a constant rate, 1/a oc 20.6 x t, and Ugr^ the semimajor axis at which 
gravitational wave radiation becomes dominant if we assume that the BBH hardens at 
the above non-hnear rate, 1/a oc 21.7 x as is actually observed in the simulations. If 
we scale our galaxy such that Tq = 37r/16 corresponds to 100 pc, and M — 1 corresponds 
to a total mass in stars in the bulge of IO^Mq, then we obtain agr-^/ugr^ — 0.3. Thus, the 
gravitational radiation stage sets in at a larger semimajor axis if hardening occurs more 
slowly than linearly. This arises from the different behaviors of the hardening timescales: 
when 1/a (x t, = \a/a\ oc 1/a, whereas when 1/a oc t^'^, oc 1/a^. However, the times 
taken to reach and approximately 1 Gyr and 70 Gyr respectively, since in the 

latter case the BBH hardens only as t^'^ rather than oc t, as in the former case. Therefore, 
despite the fact that the gravitational radiation stage is reached at a larger semimajor axis 
if hardening occurs slower than linearly, the time taken to reach this stage is longer in the 
former case. Thus, our results in this section hold for hght BHs which, as mentioned before, 
is the case with most galaxies. 

6. Summary 

In summary, our primary results are: 

• We have developed a consistent description for the behavior of the semi-major axis of 
a hard BBH. As long as the BBH mass is small compared to the total mass in stars: 
1/a oc f^, with d 1. 

• Our simplified theory gives the value of the proportionality constant of the above 
hnear relation, s — 20.6, independent of the mass of the BBH or the number of 
stars. We find, however, from numerical experiments that s depends on the number 
of stars - as the number of stars increases, s drops, but eventually converges to a 
particular value for N ^ 200, 000 stars. A possible explanation is that this reflects 
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an equilibrium between two competing processes: as the number of stars increases, 
the wandering of the binary becomes smaller, causing a larger drop in central density; 
this, however, increases the wandering of the BBH because of the reduction of the 
central restoring force on it. 

• We also find that s depends on the mass of the BBH: in our simulations for equal-mass 
BHs, the convergent value of s increases from ~ 8 to ~ 15 as the mass of each BH 
drops from 0.005 to 0.00125. We interpret this as meaning that the theoretical value 
20.6 is reached when no change in the stellar distribution function is caused by the 
BBH - this will be true only when the mass of the BBH is very small compared to the 
total mass in stars. If we could run simulations with even lower mass BBHs, we may 
approach 20.6. However, it should be remembered that we have assumed a somewhat 
uncertain form for the energy removed from the binary by each stellar interaction. 
We also note that for calculating the time at which gravitational radiation starts to 
dominate the energy loss, we do not need to know the precise value of s, since the 
fifth root of this number needs to be taken, as shown in §5. There are two difficulties 
with running simulations with very light BBHs: first, the black holes need to be much 
more massive than the stars; second, the dynamical friction timescale for light black 
holes, and consequently the time for them to become hard, is very long. 

• For light BBHs, strong disruption of the stellar density does not occur: the central 
density of the cluster does decline with time but gently. This is where we expect 
our theory to be applicable. There is no sign that a "hole" is forming at the 
center, and thus no sign that hardening would stop after the loss cone is evacuated. 
Wandering allows the BBH to interact with regions of stellar phase space it would 
not otherwise encounter. According to our experiments, if we keep the BBH center of 
mass artificially clamped to the center, hardening quickly stops after the BBH has 
interacted with the stars that could pass within roughly a distance a from it. 

• For heavy BBHs, there is stronger disruption of the stellar density and distribution 
function, and strong deviations from the linear behavior of the inverse semi-major 
axis: 1/a oc t^, with 0.5 ^ d ^ 1. This is the regime in which our theory is not 
applicable, though it gives a qualitative explanation for why d should be bracketed by 
0.5 and 1. 

• We have shown that the hardening of the BBH does not stop even when its typical 
wandering radius is smaller than the value of the semi-major axis at which it becomes 
hard, at least in the case when the BBH is fairly massive; however, the wandering 
radius is always larger than the semi-major axis at which hardening would halt 
because of loss cone depletion if the binary did not wander. The reason involves the 
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increased wandering of the binary's center of mass as the central density of the stellar 
system falls, resulting in a reduction of the restoring force keeping the binary at the 
center. To verify this in the case of light binaries would require simulations with many 
millions of particles. 

• We have calculated an upper limit on the time required for the BBHs to reach the 
stage at which energy losses are dominated by gravitational radiation emission. For 
reasonable characteristics of the stellar system, the gravitational radiation stage is 
reached within a Hubble time. In order to reach this stage, the BBH must interact 
with a total mass in stars that is of order its own mass. 

Our model uses the spherically symmetric Plummer density profile, which has a fiat 
density structure at the center. However, observations have shown that most elliptical 
galaxies are characterized by weak or strong cusps (e.g., Faber et al. 1997, Byun et al. 
1996, Lauer et al. 1995). Is it then reasonable to use the Plummer model in calculations of 
BBH coalescence? 

The Plummer model is characterized by a simple distribution function, which makes 
it possible to obtain analytic results for the phenomena discussed in this paper. Here, we 
are mainly interested in understanding the effect of the binary's Brownian motion on its 
hardening, rather than in calculating precisely the time taken for a binary in a real galaxy 
to coalesce. Despite the fact that our model for galactic bulges is unrealistic, we believe 
that our results capture generic dynamical processes that take place in more complicated 
situations for a number of reasons. First, recent simulations have shown that the merger 
of two galaxies with steep density cusps produce a galaxy with a shallow density cusp 
(Milosavljevic & Merritt 2001), owing to the transfer of energy from the binary to the 
stars; moreover, the destruction of the steep cusp takes place very rapidly (on a timescale 
comparable to the local dynamical time), before the binary becomes significantly hard 
(similar results have been obtained by Nakano and Makino 1999 a, b). Our results apply 
to the stage when the binary has become hard. Secondly, the weak cusp of the merged 
galaxy is carried by the BBH as it wanders around; whereas the stellar density and velocity 
dispersion in our calculations refer to their values in the core of the galaxy but outside the 
Keplerian rise around the black holes, since most of the encounters that remove energy from 
the BBH are with stars from outside the cusp. Thirdly, the presence of the density cusp 
would only increase the supply of stars to the binary, resulting in more efficient hardening; 
if a binary coalesces within a Hubble time in our model, it would likely coalesce within a 
shorter time in the presence of a cusp. 
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Table 1: Values of the wandering radius, Tcm, of the binary's center of mass according to our 
simplified theory (equation [13]) and numerical experiments, and the corresponding values 

of tth- 
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Table 2: Values of the wandering radius, rem, of the binary according to our simplified theory 
(equation [13]) and numerical experiments, and the corresponding values of and aic- 
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Fig. 1. — Results from a simulation of a binary black hole system with rrii — m2 — 0.00125 

and 400,000 stars. Shown arc the black hole separation, R, the reciprocal of the semi-major 
axis of the binary's orbit, a, its eccentricity, e, and the density at the center of the potential, 
p. The binary becomes hard atl/a~l/a/i = 905; it subsequently hardens at approximately 
a constant rate, s = ^(^) — 15.5. 
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Fig. 2. — The same as Figure 1, except that mi — m2 — 0.02 and N — 200, 000. The binary 
becomes hard at 1/a ~ l/o/i = 57; it subsequently hardens at a rate given approximately 
by 1/a oc t*^, with d ~ 0.5. If the center of mass of the binary is held fixed at the origin, 
hardening stops as 1/a ~ l/ojc = 130, as indicated by the lower curve of the upper right 
panel. 




Fig. 3. — The radial distance of the binary's center of mass from the origin, and the density 
at the center of the potential, for a simulation with mi — 1712 = 0.04 and N = 100, 000. The 
wandering of the BBH's center of mass increases as the central density drops. 
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Fig. 4. — (a) The line indicates those values of mi and galaxy central density, po, for which 
the BBH reaches the gravitational wave stage (a = agr) in 10 Gyr; points above this line 
correspond to combinations for which this stage is reached in less time. The stellar velocity 
dispersion a is set to 100 km/s. (b) The hne indicates those values of mi and a for which 

the BBH reaches a — agr in 10 Gyr; points above this line correspond to combinations for 
which this stage is reached in less time, po = IOOOM0PC-2. (c) The total mass in stars the 
binary must interact with until a = agr-; a = 100 km/s and po = 1000 MqPC^^. (d) dgr as a 
function of mi, a = 100 km/s and po = 1000 MqPc~^. The solid lines are for mi = 1712, the 
dashed lines for mi = 10m2. 



